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1 . Introduction 


In the analysis of spectral methods, Neumann or third type boundary value 
problems for second-order elliptic operators have received less attention than 
Dirichlet boundary value problems. The eigenvalues of a family of Chebyshev 
collocation operators related to non-Dirichlet boundary conditions were 
analyzed in [8] , while the properties of stability and convergence of such 
schemes were investigated in [3] using a general variational principle. In 
both cases, the boundary conditions are satisfied exactly by the spectral 
solution, while the differential equation is collocated at the interior nodes. 

An alternative method of imposing the boundary conditons within a 
pseudospectral scheme consists of modifying the boundary values of the first 
derivative according to the Neumann or third type conditions, during the 
evaluation of the differential operator. The equation is now collocated at 
the boundary points, too. In this way all the grid-points are treated at the 
same way by the iterative or time advancing algorithm of solution. We call 
this method the implicit treatment of the boundary conditions. 

In this paper we prove the stability and convergence of both a Legendre 
and a Chebyshev collocation scheme in which the Neumann boundary conditions 
are treated implicitly. Global error estimates are derived. Moreoever, it is 
proved that the spectral solution satisfies the boundary conditions up to an 
error which decays spectrally. Thus the spectral accuracy of the method is 
not wasted. 

Since the spectral collocation approximations of second order boundary 
value problems are usually solved by iterative techniques (see, e.g., [13]), 
we carried out an experimental analysis on the eigenvalues of the 
corresponding operators, in which the boundary conditions are imposed either 



explicitly or implicitly. The results of this investigation are also reported 
in this paper. It was found that for both the boundary treatments the 
eigenvalues are real and positive. The matrix arising from an explicit 
treatment can be preconditioned in a more natural way. However, it is shown 
in Section 5 how to build-up an effective preconditioner also in the case of 
implicitly-treated boundary conditions. 

Notations: The following notations will be used throughout the paper. 

IP N : the space of the algebraic polynomials of degree up to N in the 

variable x; 

w(x) = (1 - x 2 ) ^ : (the Chebyshev weight), 

or 

w(x) = 1: (the Legendre weight); 

2 

L (-1,1): the Hilbert space of the (classes of) Lebesgue integrable 

w 

functions v such that the norm 

* wt 0,w = (/ v 2 (x)w(x)dx) 1/2 

is finite; 

H™(-1,1): the Sobolev space of the functions v e L^(-l,l) such that their 

2 

distributional derivatives of order up to m are in L w (-l,l), 
with norm 


«vl m - ( ! / 1 [v (k) (x)] 2 w(x)dx) 1/2 J 

m,w k=0 -1 
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H^q(- 1,1): the subspace of H^(-l,l) of the functions vanishing at the 

boundary points x = ±1; 

the Hilbert space of the (classes of ) functions 

v = v(x,t) such that for almost every t e (0,T) the 

function v(«,t) e H m (-l,l) and the norm 

w 

(/ w dt ) 1/2 

0 ’ 

is finite. 

When w(x) = 1, we will drop the subscript w in all the previous notations. 


L 2 (0,T; tfVl.l)): 

w 


2. The Treatment of the Boundary Conditions 

We shall base our discussion about different treatments of the boundary 
conditions upon the following model problems: 


( 2 . 1 ) 


-u + u = f(x), 

XX * 


-1 < X < 1 


and 

( 2 . 2 ) 


( u t - u^ = f(x,t), -1 < X < 1, t > 0 

^u(x,0) » Uq(x) , -1 < x < 1 


In both cases, the solution u is assumed to satisfy the homogeneous 


boundary conditions: 
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(2.3) 


B . u = f5u + au = 0 at x = 1 
+ x 


B u = 6u + yu = 0 at x = -1 , 
x * 


for real constants a,S,Y and 6 such that 


(2.4) 


a 2 + 0 2 * 0 a$ > 0 


Y 2 + 5 2 * 0 y5 < 0. 


Under this assumption, one has by partial integration 


1 


1 


- / u udx > / (u ) dx. 
-1 XX “ _i x 


Hence, the energy method (see, e.g., [11], Vols. I and II) assures that for 

2 2 
all f e L (-1,1) there exists a unique variational solution u £ H (-1,1) 


of the boundary value problem (2.1), (2.3), such that 


flu < Hf 1 Q . 


Similarly, for all u^ e L^(-l,l) and all f e L^(0,T; L^(-l,l)) with 

2 1 

T > 0, there exists a unique variational solution u e L (0,T; H (-1,1)) O 

oo 9 

L (0,T; L (— 1 , 1) ) for the initial-boundary value problem (2.2), (2.3), such 
that, for all t < T 

»u(t)ll 2 + J Iu(t)« 2 dt < lu 0 Hq + exp(t) J If(x)l 2 dx . 

Moreover, the regularity of the solution (in the Sobolev scale H m (-l,l) or 
H™(-1*1)) increases with the regularity of the data. 
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REMARK 2.1: Weaker assumptions than (2.4) assure the well-posedness of 
the boundary value problems (2.1)-(2.3) (see, for instance, [8], Theorem 
2.1). However, we are not interested here in the minimality of our 
hypothesis, since we want to focus on the essential aspects of the treatment 
of the boundary conditions. For the same reason, we confine ourselves to very 
simple model problems, although the methods we discuss apply to general 
boundary value problems as well. 


We want to discretize (in space) equations (2.1) and (2.2) by a 
pseudospectral collocation method of Chebyshev or Legendre type. To this end, 
we look for an approximate solution u^ which is a global algebraic 

polynomial of degree N in the domain (-1,1). Moreover, we consider the 
N + 1 nodes 


(2.5) 


“1 = x. T < 


< *0 * 1 


of the Gauss-Lobatto integration rule for the Chebyshev weight 

2 — V? _ 

w(x) = (1 - x ) i or for the Legendre weight w(x) = 1 in (-1,1), (see 
[6]). If Wj , j = 0,»»»,N, are the corresponding positive weights, one has the 
identity 


( 2 . 6 ) 


i n 

/ f(x)w(x)dx = J f(x. )w. for all f e I^N-l' 

-1 1=0 2 2 


The points w ^ , j = 0,»»*,N, are the relative extrema in [-1,1] of the N-th 
Chebyshev polynomial of first kind or of the N-th Legendre polynomial. 


Since a polynomial of degree N is uniquely defined through its values at 
the nodes (2.5) we shall identify throughout the paper a polynomial of 
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degree N with the set of its values at the same nodes. Thus, if L is a 
matrix of order N + 1 and u e Lu will denote the product of the 

rp 

matrix L by the vector (u(x^) | j = O,***,!!} 1 , i.e., 


(2.7) 


Lu — i L 


u( V 


U( V 


Vue 


Given a continuous function v in [-1,1], we denote by IjjV the unique 
polynomial of degree N, interpolating v at the nodes (2.5), i.e.. 


(2.8) I N v e %, (I N v)( Xj ) = v(Xj), j = 0,»««,N. 


Some approximation properties of the operator 1^ in the Sobolev scale 
H™, m >_ 0, have been analyzed in [1] and will be used hereafter. In 

particular, there exists a constant C > 0 independent of N such that if 
the Chebyshev points are used, one has 


(2.9) ,v " I N v" 0 , w < CN" 11 «vl m>w , V u e H“(-T,T), m > V 2 , 


while if the Legendre points are used one has 


(2.9)' Iv - I N v« 0 < CN 1/2 “ m «v« m , Vue H m (— 1 , 1) , m > V 2 • 


Finally, we recall for future reference that the semi-norm 


,V V»’ { 3 l 0 " <X J )W 3 ) 



is uniformly equivalent to the norm RvRq w over (see [1], Sections 

3.1, 3.2), i.e., there exist constants > 0, C£ > 0 independent of N 

such that 

(2.1°) Cjlv.^ < .v. N>u < C 2 lv. 0j „, * u e %. 


If the boundary conditions (2.3) are of Dirichlet type (i.e., 0=6=0), 
the typical spectral collocation scheme consists of collocating the 
differential equation at the interior nodes (2.5) and setting to zero the 
solution u N on the boundary. This procedure, which we shall call the 
explicit imposition of the boundary conditions, is not restricted to Dirichlet 
boundary conditions. Thus, the boundary value problem (2.1), (2.3) can be 
approximated as 

N ^ 

U £ IEjj 

(" U xx + uN ) (x j J = f(x j ) j = 1 > * * * >N-1 

( B + u N )(x q ) = (B_ u N )(x n ) = 0, 



while the initial-boundary value problem (2.2), (2.3) can be discretized in 
space for all t > 0 as 


( 2 . 12 ) 


f u N (t) e p n 

^ U t “ u xx^ x j’ fc ^ = f(x j’ fc > j = !»•••> N-l 

< ( B + u N )(x Q ,t) = (B_ u N )(x N ,t) = 0 
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with the initial condition u N (0) = IflUg* 

For the Chebyshev collocation points, the convergence of the scheme (2.11) 

has been proven in [3], where error estimates have also been given. 

Furthermore, in [8] it has been established that the eigenvalues obtained by 

N 

replacing f with Xu in (2.11) are all real, positive and distinct. 

The unknowns to be solved for in both (2.11) and (2.12) are the values 
of u N at the interior collocation nodes and at the boundary nodes where a 
non-Dirichlet boundary condition is imposed. The algebraic system (2.11) can 
be efficiently solved by an iterative method, applied after preconditioning 
the spectral system (see, e.g., [13], [14]). We shall base our discussion of 
the computational aspects of the boundary treatment on the Richardson method, 
which is briefly recalled in Section 5. 

The differential system (2.12) can be solved by an implicit or an explicit 
time-marching method. In the first case, one has to solve at each time step a 
discrete Helmholtz equation similar to (2.11), for which one can apply one of 
the iterative procedures proposed for spectral methods. If an explicit scheme 
is used instead, the solution is advanced only at the interior collocation 
nodes. The boundary values of u^ at the new time level are subsequently 
determined in order to satisfy the boundary conditions exactly. Such values 
are obtained by solving the 2x2 system 

( <“ + M 00 )u 0 + 6d 0N u “ ‘ ’ 6 d 0j u j 
( (y + Sd NN )u N + Sd N0 u 0 ‘ - S Jj d NJ U J 
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where {d^ } is the matrix representing the spectral derivative at the 
collocation points (2.5) (an explicit formula for {d„} can be found in 
[7]). 

A completely different strategy can be followed in the process of imposing 
the boundary conditions: these are taken into account during the spectral 
evolution of the differential operator, by modifying accordingly the boundary 
values of the first derivative of the spectral solution. Precisely, assume 
that both the boundary conditions are of non-Dirichlet type, so that one can 
set $=6=1 in (2. 3) with no loss of generality. For any v e define 
the polynomial I^j(Bv x ) as follows 



Thus Ij^(Bv x ) coincides with v x at the interior nodes, but it modifies the 
boundary values of v x according to (2.3). We consider the following 
approximation of the boundary value problem (2.1), (2.3): 
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Similarly, we discretize the initial-boundary value problem (2.2), (2.3) as 
follows: 

! u N (t) e 1P N 

u^(Xj ,t) - [i n (Bu^)] x (Xj , t) = f(Xj,t), j = (),•••, N 

uN(0) = Vo* 

Note that the differential equation is now collocated at the boundary points 
also. On the other hand, the solution is not required to satisfy - and 
generally it will not satisfy - the boundary conditions exactly. However, it 
will be proved in the next sections that the boundary conditions are satisfied 
up to an error which decays spectrally with N. 

The procedure now described, first proposed by D. Gottlieb for time 
dependent problems, will be called the implicit treatment of the boundary 
conditions. All the iterative or the time-advancing methods proposed for 
solving the approximation schemes (2.11) or (2'. 12) respectively, can be 
applied in computing the solution of (2.14) or (2.15) as well. The 
computational advantage arising from an implicit treatment of the boundary 
conditions is that the iterative process of solution acts on all the grid- 
points in the domain simultaneously. Any distinction between boundary and 
interior points is avoided. 

More complex boundary conditions than (2.3), involving integro- 

differential or non-linear boundary operators, can be easily implemented in an 
implicit way, too. For instance, in [5] a far-field radiation condition of 
the type 
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3u 

3r 


= K*u , 


where K*u is a convolution operator on the far-field boundary, was 
successfully taken into account implicitly within a Fourier-Chebyshev 
collocation method for an exterior elliptic problem. 


3. Theoretical Results for the Legendre Method 

Throughout this section we assume that the collocation points (2.5) are 
the quadrature nodes of a Gass— Lobatto formula for the Legendre weight 
w(x) =1. The corresponding weights are given by 


(3.1) 


' s= , — 

j N(N + l)[L N ( Xj )J 2 


(see, e.g., [6]), where L N (x) denotes the N-th Legendre polynomial such 

that Lj^( 1 ) — 1 . 

We shall carry out an analysis for the implicit treatment of the boundary 
conditions in the case of Neumann boundary conditions, i.e., we choose 
6=6=1 and a = y = 0 in (2.3). 

The first results concern the stability and convergence properties of the 
method. 

THEOREM 3*1: Let u^ be the solution of (2.12). The following estimate 

holds: 


/ [u N (x)J 2 dx + 2 l [uJJ(x )] 2 w < / [l f(x)] 2 dx. 
-1 1=1 x J 3 -1 N 


(3.2) 
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PROOF: Equation (2.12) can be equivalently written as 


(3.3) 


-[y Bu x)J x + uN “ V’ -1 < x < 1, 


since both sides are polynomials of degree N which match at N + 1 distinct 
points. Multiply by u N and integrate over (-1,1). Since I n (Bu^) 
vanishes at the end points (see (2.11)), we have by partial integration 


(3.4) 


J I (Bu^)u^ dx + / [u^(x) d x = / I f(x)u^(x)dx. 

j Dl X X ^ ^ IN 


On the other hand, by (2*6) 


£ ** - % w j - 


whence (3.3) follows by applying the Cauchy-Schwarz inequality to the right- 
hand side of (3.4). H 


THEOREM 3.2: Let u be the solution of the boundary value problem (2.1), 
(2.3), and u N the solution of (2.12). Assume the u e H m (-l,l) with m > 
5/2. There exist two constants C " > 0, C" > 0 independent of N, u and 
u N such that 


(3.5) 


- <\ + 1"! Is - 


9 — m 9 

< <TN Bull + C~N II f I 0 
— m m-2 
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PROOF: Set V = {v e H 2 (-l,l) | v x (±l) = o} and 

= {v e | v x (±l) = O}. Let us define first a projection operator 

v v — v 

in the following way. Given v e V, denote by w^ the orthogonal projection 
of v x upon %-i^Hq(- 1 ,1) in the inner product of H*(-l,l). According 
to [12], Theorem 1.6, we have 


Iv - w^l, < CN^ + ^ m »v» , k = 0,1. 
x k — m 


If we set 


(Rjjv)(x) = y + / w(s)ds , 


where y is such that 


1 1 
/ Rjj v(x)dx = / v(x)dx, 


-1 


-1 


it is not difficult to check by a duality argument that 


(3.6) 


„k-m 


Iv - B^vl k < CN “Ivl , k = 0,1,2. 


By definition, u = R^u satisfies the equation 

+ I - + u xx ) + (“ - “) + f - 
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W /v M 

It follows from Theorem 3.1 that the difference e = u - u satisfies the 
inequality 


/ [e N (x)J 2 dx + 2 l [e"(x )] 2 w < / [f - I f J 2 dx 
-1 j=l x J J -1 N 

+ /) Kx " "xxJ 2 dx + / [« - uj 2 dx. 

According to (3.6) and (2.9) one has 
(3.7) 

/ [e N (x)] 2 dx + 2 l [e N (x )J 2 w. _< C'N 4_2m Hull 2 + C"N 5 " 2m 
—1 j=l X J J m 

On the other hand, we have by (2*10) 


If II 


2 

m-2 * 


j, K - K - < C,I n( u x - "J'o 

< c{«u - u I 2 + 1(1 - I M )(u - u )l 2 } 

— 1 x x 0 N v x x' 0 J 

< C{lu - u J 2 + CN -1 »u - u » 2 } 

— 1 x x 0 x x 1 J 

< CN 3_2m lul 2 . 


Then (3.5) follows from (3.6)-(3.7), using the triangle inequality for 

N , . N — 

u-u =(u-u)+e. 
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REHARK 3.1 : Theorem 3.2 can also be proved by a different proof, similar 

to the one which will be given in the next section for establishing Theorem 


4.2. 


The previous theorems guarantee the stability and convergence of the 
approximation (2.14) in the norm 

1 9 N-l lf 

(3.8) IvH* = (/ V (x)dx + l v” (x . )w . ) ' 2 . 

-1 j=l X J J 

Therefore, we are led to investigate the relationship between this norm and 
the usual energy norm 

Ivl. = (J v 2 (x)dx + / v 2 (x)dx) ^ 2 . 

-1 -1 x 

The two norms are clearly equivalent for polynomials of degree up to N, in 
the sense that 

(3.9) Ivl* < Ivlj < C(N)lv«*, ¥ v e %, 

where C(N) is a function of N. However, the two norms are not uniformly 

equivalent, i.e., C(N) cannot be bounded independently of N. For instance, 

2 2 

take v = L«j, the N-th Legendre polynomial. One has IL. T I + = IL„l n = 
vt N * N U 

= 2/(1 + 2N), but = 2 + 2/(1 + 2N). The numerical evaluation of 

C(N) shows that 

(3.10) C(N) » N 3 ^ 2 as N 09 


T 
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Table 1. The constant C(N) 


N 

C(N) 

C(N)/N 3/2 

2 

4. 

1.42 

4 

9.6 

1.19 

6 

16.6 

1.13 

8 

24.8 

1.10 

10 

34.0 

1.08 

12 

44.2 

1.06 

14 

55.2 

1.05 

16 

67.0 

1.05 

18 

79.6 

1.04 

20 

93.0 

1.04 


The asymptotic behavior of C(N) observed experimentally can be 
mathematically proved as follows. 

THEOREM 3.3: Let C(N) be defined by (3.9). Then 

(3.11) C(N) < CN 3/2 

PROOF : By (2.6), 

llvlj .= HvB^ - v^(l)w 0 - v^(-l)w N V v e %, 
where wq, w n are given in (3.1). Then the theorem follows immediately from 


the next lemma 
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LEMMA 3.1: There exist two positive constrants and C£ such that 

for all v e 1„ 

N-l i 

(3.12) |v x (±l)| < Cj N 5/2 lv « 0 + C 2 N 3/2 [v x ( Xj )J 2 Wj ) /2 

PROOF : Recalling the expansion of v x in terms of Legendre polynomials 

(see, e.g., [9], Appendix) one has 

1 

/ v x (x)L N (x)dx =0 

1 1 
/ v x (x)L N _ 1 (x)dx = (2N + 1) / v(x)L N (x)dx 


or equivalently 


(3.13) v x (l)w 0 + (-1)” v x (-l)w N - v x (x J )I 1J (x j )» J 


(3.14) 


; (D " 0 - <- J)M V x ( - 1)W N ' -Jj v x < 1 ‘j )L N-l (x j ) "j 


+ (2N + 1) / v(x)L„(x)dx. 
-1 N 


Then (3.12) follows using (3.1) and the Cauchy-Schwarz inequality on the 
right-hand side. H 


The spectral solution u^ of problem (2.12) is not required to satisfy 
(and generally it will not satisfy exactly) the boundary conditions (2.3). 
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However, since these are taken into account in the collocation process, one 

expects u N to satisfy approximate boundary conditions very close to the 

N 

exact ones. Lemma 3.1 yields an estimate for the values of u on the 

x 

boundary. Actually, define u = R^u as in the proof of Theorem 3.2 and set 

again e^ = u - u^, so that |e N (±l)| = |u^(±l)|. Using (3.12) for v = e^ 

and (3.7), one obtains the following result. 

THEOREM 3.4: Under the hypothesis of Theorem 3.2, the following estimate 

holds: 

(3.15) |u N (±i)| < c'N 9 ^ 2 “ m iiuii + <r'N 5-m Ufa 0 . 

1 x 1 — ra m-A 


Estimate (3.15) shows in particular that the boundary conditions are 
satisfied with spectral accuracy when they are imposed implicitly in a 
collocation scheme. 

The analysis of stability and convergence for the discrete initial 
boundary value problem (2.13) can be carried out by the same technique used in 
the proofs of Theorems 3.1 and 3.2. We omit the details of the analysis and 
we report hereafter the final result. 

THEOREM 3.5: Let the solution u of the initial-Neumann boundary value 

problem (2.2)-(2.3) satisfy the following regularity assumptions: 

u e L 2 (0,T; H m (-l,l)), u t e L 2 (0,T; H m " 2 (-l,l)) 

for a fixed T > 0 and m > 5/2. If u^ denotes the solution of the 

approximation scheme (2.13), then for all t < T 
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»u(t) - u N (t)S + (J l [ u - u^j 2 (x t)w ) V 2< C'N 2_m lu n l 

0 j = 1 A * J J U Wr~Z 

< 3 - 16 ) + exp(|-){C~N 2_m [j 1 lu(x) I 2 dx + J* lu (x)» 2 dxl ^ 

0 m 0 t 

+ C'"’N 5/2 ” m (/ C If (t) n 2 dx) 1 ^}, 

0 m 1 

for suitable constants C , C" , Q,'" independent of u, and N. 


4. Theoretical Results for the Chebyshev Method 

The most popular family of collocation points is the family of the 
Chebyshev points, which we consider in this section. The nodes (2.5) are 
given by 


(4.1) 


X. = cos 
J 


±1 

N 


j = 0, • • • ,N, 


while the corresponding weights are 


•jj j = 1 » * * * ,N-1 

ir 

W N = 2N * 

Hereafter, we shall discuss some theoretical properties of the implicit 
treatment of the Neumann boundary conditions. From now on we assume that 
0=6=1 and a = y = 0 in (2.3). 



r 



First, we prove the stability and the convergence of the time-independent 
scheme (2.1 A). 

THEOREM 4.1: Let u N be the solution of (2.14). There exists a 

constant C > 0 independent of N and such that 

(A. 3) / [u N (x)j 2 w(x)dx + Y [uj(x )] 2 w < C J 1 [l f(x) J 2 w(x)dx. 

-1 j=l X 3 2 -1 

PROOF : Equation (2.1 A) can be equivalently written as 

(A. A) -I m (Bu N ) + u N = I M f, -1 < x < 1 

v N v x'x N 

since both sides of (4.4) are polynomials of degree N which match at 
N + 1 distinct points. Following an idea due to D. Gottlieb, let us 

differentiate (4.4) with respect to x. If we set U (x) ® I^[Bu x J(x), then 
U N is a polynomial of degree N which vanishes at the boundary points and 
satisfies the collocation equation 

(4. 5) -U^( Xj ) + U N (x.) = (I N f)( Xj ), j - 1, — ,N-1. 

The stability analysis for the Chebyshev collocation approximation of the 
Dirichlet boundary value problem (see [2]) yields the estimate 
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for a constant C > 0 independent of N. Using this estimate and equation 
(4.4) we obtain the further inequality 


(4.7) 

This proves (4.3). 


Ju N n_ < cm. T fn. . 

0,w — N 0,w 


THEOREM 4.2: Let u be the solution of the boundary value problem (2.1), 

(2.3) and u N the solution of (2.14). Assume that u e H m (-l,l) with m > 

w 

5/2. There exists a constant C > 0 independent of N, u and u N such that 


(4.8) 


,u ' a \,* + ( Jj ["x ‘ u xJ (x j )w j) 1/2 < ™ 2 '” 


nun 


m.w 


PROOF: The convergence analysis for the approximation (4.5) gives the 

estimate ([2]) 


N N 

Hu - U B n + lu - ITI. < C, 
xx x O.w x O.w — 1 


l 2_m Jull + C, N 2 “ If II . 

m.w 2 m— 2,w 


(4.9) 


< C. N 2_m Bui 
— 3 ra,w 


where we have used (2.1) in order to bound the norm If I „ by the norm 

m-2,w 

w * ® y equations (2.1) and (4.4) we get 


lu - u N I n < C. N 2_m 
O.w — 4 


Bull 

m.w 
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On the other hand, the equivalence of norms (2.10) and the triangle inequality 
yields the inequality 


N-l 


K - u ") 2(x j , "j 1 SK • Vo.w + ,u x - ""'o.wl- 


whence the result, using (4.9) and (2.9). 


As for the Legendre method discussed in the previous section, it is 

possible to estimate the error on the boundary conditions produced by the 

spectral solution. We have the following result. 

THEOREM 4.3: Under the hypothesis of Theorem 4.2, there exists a 

constant C > 0 independent of N and u^ such that 

(4.10) |«£(±i)| < CN 4 " m Hu« m>w 

PROOF : For any polynomial v e IPjg, one has (see, e.g., [9]) 

1 

/ v (x)T (x)w(x)dx = 0 
_ ^ x N 

1 1 
/ v x (x)T N-1 (x)w(x)dx = N / v(x)T N (x)w(x)dx, 

whence 

2 £[v x < 1 ) + (-1)" v x (-l>] - 


(4.11) 
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(4.12) 


^■[v x (l) - (-1) N v x (-l)j = v x (x j )T N _ 1 (x j )w j 


+ N / v(x)T M (x)w(x)dx. 
-1 w 


Let u e Pjg be a primitive of the ^-projection of u x upon the space 

{v e P N-1 |v(±l) = 0}. Then (4.10) follows from (4.11) and (4.12), choosing 
, N ~ 

here v = u - u and using Theorem 4.2 in order to estimate the right-hand 
sides* wm 


As far as the evolution scheme (2.15) is concerned, the following 
convergence estimate holds. 

THEOREM 4.4: Let the solution u of the initial— Neumann boundary value 

problem (2.2), (2.3) satisfy the following regularity assumptions: 


u e L 2 (0,T; h”(- 1,1)), u t e L 2 (0,T; ^(-l.l)) 


for a fixed T > 0 and m > 5/2. Moreover, let u Q e H®(-1,1). If u 
denotes the solution of the approximation (2.15), one has for all t < T 


«u(t) - 


uN(t), o,» + (Jj K ■ u ^ 2 <x j' t)w j ) 1/2 


< 0» 2 ‘"(lu n l „ + «P (!•)[/ lu(T)lfdT 


0 ra,w ' ri 2 J i-J 


*1 
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PROOF: As in the proofs of Theorems 4.1 and 4.2 we set 


U N (x,t) = I N (Bu^)(x,t) 

and we x- differentiate the equation 

U? - I m (Bu N ) «' I N f , 
t x'x N 

which is equivalent to (2.15) . Estimate (4.13) is then a consequence of the 
convergence analysis for the Chebyshev approximation to an initial— Dirichlet 
boundary value problem (see [2], Theorem 3.3). The error on the initial data 
is 

u 0,x ' uN<0) ' u 0,x - I n( Bu x ( 0> ) - u 0,x - V*)' 


O 

The L - norm of this term can be estimated as follows: 
w 


,u 0,x - ° H(0 > , 0,w < ,u 0,x ' J N u O,x'o,w + "nK.x - B(I N “oU'o.w 

(by (2.10)) < lu 0 x - I^j u 0>x l 0>w + C.u 0jX - B(Ijj u 0 ) x l H u . 


Since both u rt and B(I„ u n ) are zero at the boundary, one has 
0,x N 0 x 


N 


,u 0,x - B(I N K,x - (I N "oU <V B J 


- CII N u 0,x " u 0 ) x , 0,w 


1 C ^ u 0,x “ V^.x^O.w + “ u 0,x _ (I N “oVo.w^ 
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By the estimate (3.7) in [1] we conclude that 


l u 0 iX -U N (0). 0i „< CN 2 -" 


The remaining part of the proof is straightforward. 


5. Computational Aspects of the Methods 
5.1 The Richardson Iterations 

The Richardson method with a finite difference preconditioning (see [13]) 
is certainly the simplest and most popular iterative method for solving 
spectral systems. We shall briefly discuss the use of this method in the 
solution of the linear systems (2.11) and (2.14). Hereafter we assume again 
that 3=5=1 in (2.3). 

The system (2.11) can be written as 


(5.1) 


l e "" - F E' 


where L E is the matrix of order N + 1 defined by the relation (recall 
(2.7)): 


(-B + v)(x q ) 


(5.2) L „ v - I (-v + v)(x.) I j = !,•••, N-l, for all v e 


XX 


%> 


(B v)(x ) 

- N 


T 
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and F e = (0,f (xj) , • •• ,f (x^_j) ,0)^. An approximate inverse Ag of Lg, to 
be used as a preconditioner, can be built up by low order finite differences 
at the nodes (2.5) as follows: 


(5.3) 


\ v 



for all v e IPjg, where hj = Xj - Xj + p Thus, a one-sided finite difference 
approximation of the boundry conditon is imposed at the boundary nodes, while 
centered differences are used at the interior. 

The system (2.14) is represented as 


(5.4) 



where Lj is the matrix of order N+l defined by the relation 

(5.5) Lj v - {["V Bv x ) x + vl(x j^0<j<N for 311 V e % 

and Fj = {f(Xj ) Jo<j<N* Preconditioning this matrix is a more delicate 
matter than preconditioning the matrix (5.2). In analogy with (5.3), one 
could consider the matrix Aj defined as 
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name ly , the differential operator is discretized also at the boundary nodes by 
a centered difference formula, and the boundary conditions are used in 
eliminating the auxiliary nodes outside the domain. Such a matrix exhibits 
very poor preconditioning properties for the matrix (5.5). This can be 
explained by considering the structure of the spectrum of Lj in the case of 
Neumann boundary conditions. The eigenvalue 1 has double multiplicity, the 
corresponding eigenfunctions being v = 1 and v = <f> N (the N-th Chebyshev or 
Legendre orthogonal polynomial). Actually, the Gauss— Lobatto quadrature nodes 
are such that 4> N>x (Xj) = 0 for j = 1,***,N-1. Hence 


V B V*> 


0 . 


On the contrary, is not an eigenfunction for the matrix (5.6) and 1 is a 
simple eigenvalue. 

However, it is easy to build up a finite difference approximation of the 
operator (2.1) which for the Neumann boundary conditions has 1 as a double 
eigenvalue with eigenfunctions v = 1 and v = <k, T . At each interior node 
Xj (j-1 , • •• ,N-1) , define the differentiation formula 

— 
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(5.7) a u(x x ) + b j u ( x j) + c j u ^ x j+p “ u x ^ x j^ 

by the conditions of being first order accurate and of satisfying the identity 

VN^j-l 5 + b j W + C jV X j+l ) = °* 

If the Chebyshev points (4.1) are used, then <J>^(Xj) = T^(x^) = (-1)^, hence 
bj = 0 and (5.7) is the centered difference approximation of the first 
derivative. 

After computing the numbers a j , bj and Cj we define the matrix Aj as 
follows 



The matrix Aj is a five-diagonal approximation of the operator Lj, which 
has the required spectral properties. 

The Richardson method is applied in solving (5.1) and (5.4) in the 
following form: 


N,k+1 N,k . k .-1 _ N,k N 
u * = u +aA (F-Lu) 


(5.9) 


k > 0, 
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where a >0 is an acceleration parameter. In (5.9), A, F and L stand 

for Ajj, Fg, Lg if the system (5.1) is to be solved, and for Aj, F , Lj if 

the system (5.4) is considered instead. Under the assumption that the 

matrix A is diagonalizable, the method is convergent (with a proper 

choice of the acceleration parameter) provided all the eigenvalues of A“^L 

have strictly positive real parts. For the pure Dirichlet boundary value 

problem this is true since it has been proven ([10]) that the eigenvalues of 
—1 2 

A L lie in the interval [1 ,tt /4]. For the general boundary conditions 

(2.3), the behavior of the eigenvalues has been investigated numerically and 
will be discussed hereafter. 

The convergence of the Richardson method is crucially influenced by the 
choice of the parameter a . Several strategies have been proposed ([13], 

[14], [10]). The simplest and most effective one consists of computing a k 
by minimizing some Hilbertian norm of the residual r k+1 = F - Lu N,k+1 . 
Assume that the Hilbertian norm is defined by the inner product 


((u,v)) = l h u v , 
ij J 2 


where H = {h^} is a symmetric positive definite matrix. Then, the 
resulting expression for o is 


(5.10) 


flk _ ((r k , LA“ X r k )) # 

((LA -1 r k , LA -1 r k )) 


The iterative procedure (5 .9)— (5. 10) (called the Minimal Residual Richardson 
Method) will not break down if a k remains bounded away from 0. This occurs 
if the symmetric part of the matrix HLA” 1 is positive definite. 
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Unfortunately, the numerical computations described in the next subsection 
show that whenever the boundary conditions are not of Dirichlet type few 
eigenvalues of HLA“* + (HLA”^) T are of negative sign. 

Thus one has to resort to iterative methods which converge even if the 
symmetric part of the matrix of the system is indefinite. Among them, the 
algorithm Orthodir (see [15]) seems particularly apt to be used with spectral 
methods. Setting B = A - * L, the algorithm is defined as 


u 


N,k+1 


N,k 

u 


+ a 


k 

P 


where the descent direction p is given by 


k k-1 ^7^ i 

P = -J, B k . p J , 

j=o kJ 


(B 2 p k 1 , Bp 3 ) 
(Bp j , Bp j ) 


kj 


and a is chosen by the minimal residual strategy. 

A truncated version, consisting of setting =0 if j < k-1, is 

generally preferred, although the convergence is not assured in this case. 
Since one step of Orthodir requires twice as many operations as one step of 
Richardson^s method, it is convenient to execute a few Orthodir iterations 
only when the Richardson method breaks down, then going back to the original 
method. See [5] for a successful application of this strategy. 


5.2 The Behavior of the Eigenvalues 

The eigenvalues of the spectral matrices L, HL + (HL) T and the 

preconditioned matrices A~^L, HLA~^ + (HLA~^)^ were computed by EISPACK 

k 

routines. The computation was carried out for N = 2 , k = 2, •••,6 and for 
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the boundary conditions a = y = 0 (pure Neumann) a = 1, y = 0 (mixed 
Neumann-Robin) and a = y = 5 (pure Robin). The eigenvalues were found 
scarcely sensitive to the kind of boundary conditions, the qualitative 
behavior being the same for the three conditions considered. The following 
discussion refers to the pure Neumann problem. 

a) Chebyshev explicit method ((5.1), (4.1)) 

The eigenvalues of the matrix Lg defined in (5.2) are real, positive, 
distinct and bounded from below by 1. The eigenvalues of the preconditioned 
matrix L^, with defined in (5.3), are also real positive and 
distinct. Moreover, they lie in the interval [1,tt 2 /4], as shown in Table 2. 


Table 2. Largest eigenvalues for the 
Chebyshev "explicit" method 


N 

l e 

a e* l e 

4 

24. 

1.75 

8 

231. 

2.10 

16 

3242. 

2.29 

32 

50,208. 

2.38 

64 

701,902. 

2.44 


Let H be the matrix associated with the Chebyshev discrete inner product, 
i.e., H *» diagfwQ , • • • ,w^} where the Wj's are defined in (4.2). The 
symmetric part of the matrices HL^ and HL^ is indefinite. In both 
cases two eigenvalues are negative and their largest modulus is of the order 
of the largest positive eigenvalues. Unlike A^ 1 L £ , the eigenvalues of 
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hl e + 

eigenvalue growing like O(N^) • A similar behavior occurs if H is the 

identity matrix* 

b) Chebyshev implicit method ((5*4), (4*2)) 

The eigenvalues of the matrix Lj- defined in (5*5) are real and positive, 
as it easily follows from (4.5). Moreover, they are all simple, except the 

smallest eigenvalues 1 which is double, as already pointed out. If the 

matrix (5.6) is used as a preconditioner, the largest eigenvalues of the 

matrix L^. remains bounded, but the smallest eigenvalue tends to zero 

as 0(N~3). 

If the matrix (5.8) is used instead, the smallest eigenvalue of the matrix 
A” 1 Lj. is 1 with double multiplicity, while the modulus of the largest 

eigenvalue, although not bounded independently of N, grows like 0 (N^^). 


Table 3. Largest modulus of the eigenvalues 

for the Chebyshev "implicit" method. 
The matrix Aj Is defined in (5.8). 


N 

C(N) 

a I L Lj 

4 

19. 

5. 

8 

214. 

14.6 

16 

3174. 

38. 

32 

49,938. 

107. 

64 

700,893. 

342. 


(HL E A £ ) are not bounded independently of N, the largest 
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The eigenvalues of L^ are all real and positive, except the two largest 

eigenvalues which for N ^ 32 are complex conjugate (with positive real 
parts) . 

The symmetric part of the matrices HL X and HL^ A” 1 are indefinite. 
The number of negative eigenvalues for HL^ + (HLj.) 7 is 2 in all the cases 
under investigation, while for HLj. a" 1 + (HL ]; a" 1 ) 1 this number grows slowly 
with N (it is 6 for N = 64). 

c) Legendre explicit method ((5.1), (3.1)) 

The eigenvalues of Lg and Ag^ Lg behave qualitatively as the 
eigenvalues of the corresponding matrices for the Chebyshev method. Table 4 
contains the largest eigenvalues for the two matrices, the smallest eigenvalue 
being 1 in all cases. 


Table 4. Largest eigenvalues for the Legendre 
"explicit" method 


N 

l e 

a e 1 l e 

4 

22. 

1.6 

8 

162. 

1.9 

16 

1978. 

2.17 

32 

28,639. 

2.31 

64 

432,449. 

2.41 
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For the Legendre method too, the symmetric parts of HLg and HL^ A E * (where 
H = diagfwQ,*** ,w^}, Wj being defined in (3.1)), are indefinite. 

d) Legendre implicit method ((5.4), (3.1)) 

The eigenvalues of the matrix Lj are real, positive and bounded from 
below by 1. Such properties can be easily proved using the identity (3.4). 
The same properties are shared by the eigenvalues of the preconditioned matrix 
Lj, with Aj defined in (5.8); in this case the largest eigenvalue grows 
like 0(N 3 / 2 ). 


Table 5. Largest eigenvalues for the 
Legendre "implicit" method 


N 

l e 

Ai 1 Lj 

4 

17. 

5.6 

8 

146. 

18.7 

16 

1921. 

55. 

32 

28,424. 

148. 

64 

431,676. 

384. 


The symmetric parts of the matrices HLj- 


and 


HL 


I 



are indefinite 
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